Remove["Global`*"]
data =Transpose[(({{1, 2, 3, 4, 5}, {1, 3, 2, 4, 2}}))];
deg = 4;
p4[x_] =C[0]+
	C[1](x-data[[1]][[1]])+
	C[2](x-data[[1]][[1]])(x-data[[2]][[1]])+
	C[3](x-data[[1]][[1]])(x-data[[2]][[1]])(x-data[[3]][[1]])+
	C[4](x-data[[1]][[1]])(x-data[[2]][[1]])(x-data[[3]][[1]])(x-data[[4]][[1]]);

const = Table[C[i], {i, 0,deg}];
sol=Solve[Table[p4[data[[i]][[1]]]==data[[i]][[2]],{i,1,deg+1}],const];

p[x_]=Fit[data,{1,x,x^2,x^3,x^4},x];

Print["Our Vandermond formulation: ",Simplify[p4[x]/.sol[[1]]]]
Print["Fit function:               ",p[x]]
Print["Are they the same? ",Expand[p4[x]/.sol[[1]]]==Expand[p[x]]]

pl=Plot[p4[x]/.sol[[1]],{x,data[[1]][[1]],data[[deg+1]][[1]]},PlotLegends->{"p4(x)"}];
lp=ListPlot[data,PlotStyle->{Red,PointSize[Large]},PlotLegends->{"data"}];
Show[{pl, lp}]

Print["Compare with the result given by Fit. Estimate y(4.2) and y'(4.2). Compare with y(4.2) from 2.2"]
Print["(From 2.2) Interpolation[data, 4.2]=",Interpolation[data, 4.2] ](*from 2.2*)
Print["Estimated y(4.2) with Fit. p[4.2]=",p[4.2]] (*Estimated y(4.2) *)
Print["Estimated y(4.2) with our function. p4[4.2]=",p4[4.2] /.sol[[1]]](* Estmited y(4.2) *)
p4d[x_] = D[p4[x],x]/.sol[[1]]; (* cool derivate function *)
Print["Estimated y'(4.2) with our p4d[4.2]=",p4d[4.2]]